load packages and functions

source('mytheme.R')
# model version
modelversion = 'log_rstan'
rstanmodelPath = 'modelrlt'
modelPath = paste0(rstanmodelPath, '/models/', modelversion)

Merge the model Result data

1 load data

AllExpData = read.csv(paste0("../data/AllValidData.csv"))
dur <- sort(unique(AllExpData$curDur))

AllExpData$WMSize <- factor(AllExpData$WMSize, labels = c("low", "medium",  "high")) 
# 1: 500ms,  2: 2500, 3 : 2000ms the mean reaction time of WM test 
AllExpData$gap <- factor(AllExpData$gap,  labels = c('short','long', 'not sure'))
AllExpData[which(AllExpData$Exp == 'Exp4'),"Exp"] = "Exp4a"
AllExpData[which(AllExpData$Exp == 'Exp5'),"Exp"] = "Exp4b"

2 Corrct rate

WMCrr2 <- ggplot(meanForPlot, aes(WMSize, mean_WMCrr, ymin = mean_WMCrr - se_WMCrr, ymax = mean_WMCrr + se_WMCrr,
                                  group =Exp, color = Exp, fill = Exp))+ 
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.2,  position = position_dodge(width = 0.2)) +
  coord_cartesian(ylim = c(0.5, 1)) +
  colorSet5+
  labs(x = "Memory load", y = "Mean accuracy in WM task") +
  theme_new 
WMCrr2

ggsave(paste0(getwd(), "/figures/WMCrr2.png"), WMCrr2, width = 4, height = 4)
### generate WM correct rates
AllExpData$WMCrr <- AllExpData$TPresent == AllExpData$WMRP
m_wmp<- dplyr::group_by(AllExpData, Exp, WMSize, NSub) %>% 
  dplyr::summarize(m_WMCrr = mean(WMCrr), n =n(), se_WMCrr = sd(WMCrr)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.

3 comparison between Exp.4a and Exp.4b)

3.1 correct rate

#plot WM correct rates
dplyr::group_by(AllExpData%>%filter(Exp %in% c('Exp4a', 'Exp4b')), Exp, WMSize, NSub, gap) %>% 
  dplyr::summarize(m_WMCrr = mean(WMCrr), n = n(), se_WMCrr = sd(WMCrr)/sqrt(n-1)) -> correctrate_gap
## `summarise()` has grouped output by 'Exp', 'WMSize', 'NSub'. You can override
## using the `.groups` argument.
write.csv(correctrate_gap, paste0(modelPath, '/rlt/correctrate_gap.csv'))

correctrate_gap%>%
  dplyr::group_by(Exp, WMSize, gap)%>%
  dplyr::summarize( n = n(),
             mean_WMCrr = mean(m_WMCrr), se_WMCrr = sd(m_WMCrr)/sqrt(n-1) ) -> meanForPlot_gap
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
WMCrr_gap <- ggplot(meanForPlot_gap, aes(WMSize, mean_WMCrr, ymin = mean_WMCrr - se_WMCrr, ymax = mean_WMCrr + se_WMCrr, group =interaction(Exp,gap), color = gap, shape = Exp, linetype = Exp))+ 
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.2,  position = position_dodge(width = 0.2)) +
  coord_cartesian(ylim = c(0.5, 1)) +
  colorSet5+
  labs(x = "Memory load", y = "Mean accuracy in WM task") +theme_new 
WMCrr_gap

ggsave(paste0(getwd(), "/figures/WMCrr_gap.png"), WMCrr_gap, width = 4, height = 4)

3.2 observed RP

AllExpData%>% filter(Exp %in% c('Exp4a', 'Exp4b')) %>%
  dplyr::group_by(Exp, curDur, WMSize, NSub, gap) %>% 
  dplyr::summarize(n = n(),
                   m_repDur = mean(repDur),
                   se_repDur = sd(repDur)/sqrt(n-1)) ->RP_obs_gap_bias
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize', 'NSub'. You can
## override using the `.groups` argument.
RP_obs_gap_bias$rep_bias = RP_obs_gap_bias$m_repDur - RP_obs_gap_bias$curDur


ezANOVA(data = RP_obs_gap_bias%>%filter(Exp =='Exp4b'), dv= rep_bias, wid=NSub, within= .(curDur, gap, WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: "curDur" will be treated as numeric.
## Warning: You have removed one or more levels from variable "gap". Refactoring
## for ANOVA.
## Warning: There is at least one numeric within variable, therefore aov() will be
## used for computation and no assumption checks will be obtained.
## $ANOVA
##              Effect DFn DFd          F            p p<.05          ges
## 1            curDur   1  15 78.8694018 2.325297e-07     * 8.072483e-01
## 2               gap   1  15  0.3640022 5.553058e-01       5.913116e-04
## 3            WMSize   2  30  9.8759805 5.066007e-04     * 6.980625e-02
## 4        curDur:gap   1  15 10.6242360 5.279692e-03     * 1.079097e-02
## 5     curDur:WMSize   2  30 13.1667406 7.858144e-05     * 2.343741e-02
## 6        gap:WMSize   2  30  3.1032330 5.956988e-02       2.529934e-03
## 7 curDur:gap:WMSize   2  30  0.1127064 8.937894e-01       7.606329e-05
RP_obs_gap  <- ggplot( data = RP_obs_gap_bias%>%
  dplyr::group_by(Exp, curDur, WMSize, gap) %>% 
  dplyr::summarize(m_m_repDur = mean(m_repDur),
                   m_se_repDur = mean(se_repDur)), aes(x = curDur, y = m_m_repDur,  color=as.factor(WMSize), shape = gap)) +
  geom_point()+
  geom_line()+
  geom_abline(slope=1, intercept=0)+
  facet_grid(cols = vars(Exp)) +
  labs(x="Sample intervals (s)", y="Reproduction (s)", shape=" ", color = "Memory Load")+
  theme_new+colorSet3+
   scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override
## using the `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_obs.png"), RP_obs_gap, width = 6, height = 4)

RP_obs_gap

3.3 central tendency effect (slope) and IPs (linear regression)

#Observed Indifference Point
obs_model <- function(df) {
  lm(repDur ~ curDur, data = df)
}

obs_InP <- AllExpData %>% filter(Exp %in% c('Exp4a', 'Exp4b')) %>%
  dplyr::group_by(NSub, Exp, WMSize, gap) %>% nest()  %>%
  mutate(model = map(data, obs_model)) %>%  # linear regression
  mutate(slope = map(model, broom::tidy)) %>%  # get estimates
  unnest(slope, .drop = TRUE) %>% # remove raw data
  select(-std.error,-statistic, -p.value) %>%  # remove unnessary columns
  spread(term, estimate) %>%   # spread stimates
  dplyr::rename(Intercept = `(Intercept)`, slope = curDur)  # rename columns
## Warning: The `.drop` argument of `unnest()` is deprecated as of tidyr 1.0.0.
## All list-columns are now preserved.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
obs_InP$model = NULL
obs_InP$data = NULL
obs_InP$inP = obs_InP$Intercept /(1-obs_InP$slope)

write.csv(obs_InP, paste0(modelPath, '/rlt/obs_InP_gap.csv'))
mobs_InP = obs_InP %>% group_by(Exp, WMSize, gap)%>%
  dplyr::summarise(n=n(),
                   m_Intercept = mean(Intercept),
                   se_Intercept= sd(Intercept)/sqrt(n-1),
                   m_inP = mean(inP),
                   se_inP = sd(inP)/sqrt(n-1),
                   m_slope = mean(slope),
                   se_slope = sd(slope)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
plt_InP_linear_gap<- ggplot(data = mobs_InP, aes(x=WMSize, y=m_inP, group = gap, color = gap))+
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.3,  aes(ymin = m_inP - se_inP, ymax = m_inP + se_inP), position = position_dodge(width = 0.2)) +theme_new+
  labs(colour = "Gap")+colorSet3+
  facet_wrap(~Exp)+
  xlab(' ')+ylab("indifference point (s)")+guides(shape="none")+
  theme(legend.position = "top")

ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_InP_linear_gap.png"), plt_InP_linear_gap, width = 5, height = 5)

plt_InP_linear_gap

ezANOVA(data = obs_InP%>%filter(Exp =='Exp4b'), dv= inP, wid=NSub, within= .(gap, WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: You have removed one or more levels from variable "gap". Refactoring
## for ANOVA.
## $ANOVA
##       Effect DFn DFd        F           p p<.05         ges
## 2        gap   1  15 2.942083 0.106876928       0.006249289
## 3     WMSize   2  30 7.307023 0.002598632     * 0.069252661
## 4 gap:WMSize   2  30 1.289560 0.290220752       0.006533781
## 
## $`Mauchly's Test for Sphericity`
##       Effect         W           p p<.05
## 3     WMSize 0.3727819 0.001000422     *
## 4 gap:WMSize 0.5400426 0.013396643     *
## 
## $`Sphericity Corrections`
##       Effect       GGe      p[GG] p[GG]<.05       HFe       p[HF] p[HF]<.05
## 3     WMSize 0.6145458 0.01069076         * 0.6414052 0.009678944         *
## 4 gap:WMSize 0.6849515 0.28394478           0.7306789 0.285474638
plt_RP_slope_linear_gap<- ggplot(data = mobs_InP, aes(x= WMSize, y=m_slope, group = gap,color = gap, shape = gap))+
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.3,  aes(ymin = m_slope - se_slope, ymax = m_slope + se_slope), position = position_dodge(width = 0.2)) +theme_new+
  facet_wrap(~Exp)+
  labs(colour = "Gap")+colorSet3+
  xlab(' ')+ylab("Slope of reproduction (s)")+
  theme(legend.position = "top")

ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_RP_slope_linear_gap.png"), plt_RP_slope_linear_gap, width = 5, height = 5)

plt_RP_slope_linear_gap

ezANOVA(data = obs_InP%>%filter(Exp =='Exp4b'), dv= slope, wid=NSub, within= .(gap, WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: You have removed one or more levels from variable "gap". Refactoring
## for ANOVA.
## $ANOVA
##       Effect DFn DFd          F            p p<.05          ges
## 2        gap   1  15 10.6363088 5.259698e-03     * 1.264088e-02
## 3     WMSize   2  30 13.1535570 7.913522e-05     * 2.751879e-02
## 4 gap:WMSize   2  30  0.1142223 8.924457e-01       9.048962e-05
## 
## $`Mauchly's Test for Sphericity`
##       Effect         W         p p<.05
## 3     WMSize 0.8539487 0.3311481      
## 4 gap:WMSize 0.9665181 0.7878977      
## 
## $`Sphericity Corrections`
##       Effect       GGe       p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 3     WMSize 0.8725613 0.000189114         * 0.9778272 9.205975e-05         *
## 4 gap:WMSize 0.9676028 0.886521995           1.1084480 8.924457e-01

4 load model results

4.1 check model parameters

### Average Parameters 
mm_Baypar <- dplyr::group_by(Bayparlist, Exp) %>%
  dplyr::summarize( m_sig_s2 = mean(sig_s2),
                    m_sig_pr2_log = mean(sig_pr2_log),
                    m_ks= mean(ks), 
                    m_kr = mean(kr), 
                    m_ls = mean(ls), 
                    m_mu_pr_log= mean(mu_pr_log),
                    n= n(),
                    se_sig_s2 = sd(sig_s2)/sqrt(n-1),
                    se_sig_mn2 = sd(sig_mn2)/sqrt(n-1),
                    se_sig_pr2_log = sd(sig_pr2_log)/sqrt(n-1),
                    se_ks = sd(ks)/sqrt(n-1),
                    se_kr = sd(kr)/sqrt(n-1),
                    se_ls =sd(ls)/sqrt(n-1),
                    se_mu_pr_log = sd(mu_pr_log)/sqrt(n-1)) 
mm_Baypar

4.2 Average Parameters

mBaypar_sub <- dplyr::group_by(Bayparlist, Exp, NSub) %>%
  dplyr::summarize( m_sig_s2 = mean(sig_s2),
                    m_sig_pr2_log = mean(sig_pr2_log),
                    m_ks= mean(ks), 
                    m_kr = mean(kr), 
                    m_ls = mean(ls), 
                    m_mu_pr_log= mean(mu_pr_log),
                    n= n(),
                    se_sig_s2 = sd(sig_s2)/sqrt(n-1),
                    se_sig_mn2 = sd(sig_mn2)/sqrt(n-1),
                    se_sig_pr2_log = sd(sig_pr2_log)/sqrt(n-1),
                    se_ks = sd(ks)/sqrt(n-1),
                    se_kr = sd(kr)/sqrt(n-1),
                    se_ls =sd(ls)/sqrt(n-1),
                    se_mu_pr_log = sd(mu_pr_log)/sqrt(n-1)) 
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
fig_ks_subj = ggplot(mBaypar_sub, aes(Exp, m_ks, color = Exp)) +
  geom_boxplot(position = position_dodge()) +
  geom_jitter(shape=16, position=position_jitter(0.2))+
  #facet_wrap(~Exp)+
  theme_new+ colorSet5 +
  theme(strip.background = element_blank()) + # remove subtitle background
  labs(x = ' ', y = 'scale factor Ks')+ 
  theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_ks_subj.png"), fig_ks_subj, width = 6, height = 6)
fig_ks_subj

fig_kr_subj = ggplot(mBaypar_sub, aes(Exp, m_kr, color = Exp)) +
  geom_boxplot(position = position_dodge()) +
  geom_jitter(shape=16, position=position_jitter(0.2))+
  #facet_wrap(~Exp)+
  theme_new+ colorSet5 +
  theme(strip.background = element_blank()) + # remove subtitle background
  labs(x = ' ', y = 'scale factor Kr')+ 
  theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_kr_subj.png"), fig_kr_subj, width = 6, height = 6)
fig_kr_subj

fig_ls_subj = ggplot(mBaypar_sub, aes(Exp, m_ls, color = Exp)) +
  geom_boxplot(position = position_dodge()) +
  geom_jitter(shape=16, position=position_jitter(0.2))+
  #facet_wrap(~Exp)+
  theme_new+ colorSet5 +
  theme(strip.background = element_blank()) + # remove subtitle background
  labs(x = ' ', y = 'scale factor ls')+ 
  theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/mfig_ls_subj.png"), fig_ls_subj, width = 6, height = 6)
fig_ls_subj

5 Prediction results

#load prediction
AllDat_predY <- read.csv(paste0(modelPath, "/rlt/AllDat_predY_",modelversion,".csv"))
AllDat_predY$WMSize <- as.factor(AllDat_predY$WMSize)
levels(AllDat_predY$WMSize) = c("low", "medium",  "high")
AllDat_predY$pred_Bias = AllDat_predY$mu_r - AllDat_predY$curDur
AllDat_predY$predErr = AllDat_predY$mu_r - AllDat_predY$repDur
AllDat_predY$relatErr = AllDat_predY$predErr / AllDat_predY$repDur
AllDat_predY[which(AllDat_predY$Exp == "Exp4"),"Exp"] = "Exp4a"
AllDat_predY[which(AllDat_predY$Exp == "Exp5"),"Exp"] = "Exp4b"
#calculate the mean reproduction biases for the five given intervals for all subjects
mpredY_sub <- dplyr::group_by(AllDat_predY, curDur, Exp, NSub, WMSize) %>%
  dplyr::summarize(n = n(), 
                   m_repDur = mean(repDur), 
                   sd_repDur = sd(repDur),
                   m_mu_r = mean(mu_r),
                   m_sig_r = mean(sig_r),
                   m_wp = mean(wp),
                   se_wp = sd(wp)/sqrt(n-1),
                   log_lik =mean(log_lik),
                   cv =sd_repDur/ m_repDur,
                   pred_cv = mean(sig_r/mu_r),
                   predRP_err = mean(m_mu_r-m_repDur),
                   predVar_err = mean(m_sig_r-sd_repDur),
                   predRP_rerr = mean(abs(m_mu_r-m_repDur)/m_repDur),
                   predVar_rerr = mean(abs(m_sig_r-sd_repDur)/sd_repDur),
                   predcv_err = pred_cv-cv,
                   predcv_rerr = mean(abs(pred_cv-cv)/cv))
## `summarise()` has grouped output by 'curDur', 'Exp', 'NSub'. You can override
## using the `.groups` argument.
mpredY_sub_split <-split(mpredY_sub %>%select(c('Exp', 'NSub', 'WMSize', 'm_repDur', 'curDur')), mpredY_sub$curDur) 
mpredY_sub_cv_split <-split(mpredY_sub %>%select(c('Exp', 'NSub', 'WMSize', 'cv', 'curDur')), mpredY_sub$curDur) 

mpredY_sub_split1 <-split(mpredY_sub %>%select(c('Exp', 'NSub', 'WMSize', 'curDur' ,'m_repDur')), mpredY_sub$WMSize) 
mpredY_sub_cv_split1 <-split(mpredY_sub %>%select(c('Exp', 'NSub', 'WMSize','curDur', 'cv')), mpredY_sub$WMSize) 
mpredY_sub_jasp_RP = NULL
mpredY_sub_jasp_cv = NULL
mpredY_sub_jasp_RP1 = NULL
mpredY_sub_jasp_cv1 = NULL
for(i in 1: length(mpredY_sub_split)){
  temp = mpredY_sub_split[[i]]
  curDur = unique(temp$curDur)
  temp$curDur = NULL
  colnames(temp) = c('Exp', 'NSub', 'WMSize', paste0('m_repDur_', curDur))
  
  temp_cv = mpredY_sub_cv_split[[i]]
  curDur = unique(temp_cv$curDur)
  temp_cv$curDur = NULL
  colnames(temp_cv) = c('Exp', 'NSub', 'WMSize', paste0('cv_', curDur))
  if(i == 1){
    mpredY_sub_jasp_RP = temp
    mpredY_sub_jasp_cv = temp_cv
  }
  else{
    mpredY_sub_jasp_RP = left_join(mpredY_sub_jasp_RP, temp, by=c("Exp",  "NSub", "WMSize"))
    mpredY_sub_jasp_cv = left_join(mpredY_sub_jasp_cv, temp_cv, by=c("Exp",  "NSub", "WMSize"))
  }
  
}
for (i in 1: length(mpredY_sub_split1)){
  temp1 = mpredY_sub_split1[[i]]
  WMSize = unique(temp1$WMSize)
  temp1$WMSize = NULL
  colnames(temp1) = c('Exp', 'NSub', 'curDur', paste0('m_repDur_', WMSize))
  
  temp_cv1 = mpredY_sub_cv_split1[[i]]
  WMSize = unique(temp_cv1$WMSize)
  temp_cv1$WMSize = NULL
  colnames(temp_cv1) = c('Exp', 'NSub', 'curDur', paste0('cv_', WMSize))
  if(i == 1){
    mpredY_sub_jasp_RP1 = temp1
    mpredY_sub_jasp_cv1 = temp_cv1
  }
  else{
    mpredY_sub_jasp_RP1 = left_join(mpredY_sub_jasp_RP1, temp1, by=c("Exp",  "NSub", "curDur"))
    mpredY_sub_jasp_cv1 = left_join(mpredY_sub_jasp_cv1, temp_cv1, by=c("Exp",  "NSub", "curDur"))
  }
}

write_csv(mpredY_sub_jasp_RP, paste0(modelPath, '/rlt/RP_Bias_jasp.csv'))
write_csv(mpredY_sub_jasp_cv, paste0(modelPath, '/rlt/cv_jasp.csv'))
write_csv(mpredY_sub_jasp_RP1, paste0(modelPath, '/rlt/RP_Bias_jasp1.csv'))
write_csv(mpredY_sub_jasp_cv1, paste0(modelPath, '/rlt/cv_jasp1.csv'))
write_csv(dplyr::group_by(mpredY_sub, curDur, NSub) %>%
  dplyr::summarize(m_cv = mean(cv))%>%spread(curDur, m_cv), paste0(modelPath, '/rlt/m_cv.csv'))
## `summarise()` has grouped output by 'curDur'. You can override using the
## `.groups` argument.
mpredY_sub$RP_bias = mpredY_sub$m_repDur -mpredY_sub$curDur
mpredY_sub_new <- dplyr::group_by(mpredY_sub, curDur, Exp, NSub) %>%
  dplyr::summarize(m_RP_bias = mean(RP_bias))%>% spread(curDur, m_RP_bias)
## `summarise()` has grouped output by 'curDur', 'Exp'. You can override using the
## `.groups` argument.
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp1'), paste0(modelPath, '/rlt/RP_Bias_exp1.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp2'), paste0(modelPath, '/rlt/RP_Bias_exp2.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp3'), paste0(modelPath, '/rlt/RP_Bias_exp3.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp4a'), paste0(modelPath, '/rlt/RP_Bias_exp4a.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp4b'), paste0(modelPath, '/rlt/RP_Bias_exp4b.csv'))


mpredY_sub_WMsize <- dplyr::group_by(mpredY_sub, WMSize, Exp, NSub) %>%
  dplyr::summarize(m_RP_bias = mean(RP_bias))%>% spread(WMSize, m_RP_bias)
## `summarise()` has grouped output by 'WMSize', 'Exp'. You can override using the
## `.groups` argument.
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp3'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp3.csv'))
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp4a'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp4a.csv'))
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp4b'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp4b.csv'))
#### predicted data
m_predY <- mpredY_sub%>% 
  dplyr::group_by(Exp, curDur, WMSize) %>% 
  dplyr::summarize(m_m_repDur = mean(m_repDur),
                   m_sd_repDur = mean(sd_repDur),
                   m_m_sig_r =mean(m_sig_r),
                   m_m_mu_r = mean(m_mu_r),
                   m_m_wp = mean(m_wp),
                   n = n(),
                   m_se_wp = sd(se_wp)/sqrt(n-1),
                   log_lik =mean(log_lik),
                   mpredRP_err = mean(predRP_err),
                   mpredVar_err = mean(predVar_err),
                   mpredRP_rerr = mean(predRP_rerr),
                   mpredVar_rerr = mean(predVar_rerr),
                   cv= mean(cv),
                   pred_cv = mean(pred_cv),
                   mpredcv_err = mean(predcv_err),
                   mpredcv_rerr = mean(predcv_rerr))
m_predY_acc =  mpredY_sub%>% 
  dplyr::group_by(Exp) %>% 
  dplyr::summarize(mpred_rerr = mean(predRP_rerr)*100,
                   mpredVar_rerr = mean(predVar_rerr)*100, 
                   mpredcv_rerr = mean(predcv_rerr)*100)
m_predY_acc

6 WAIC and LOO-CV

#extract waic and loo-cv from parameter list
m_WAIC <- dplyr::group_by(Bayparlist, Exp) %>%
  dplyr::summarize(n =n(),
                   m_looic = mean(looic),
                   m_waic = mean(waic),
                   se_waic = sd(waic)/sqrt(n-1),
                   se_looic = sd(looic)/sqrt(n-1),
                   m_p_loo = mean(p_loo),
                   m_elpd_loo = mean(elpd_loo),
                   m_se_looic = mean(se_looic),
                   m_se_p_loo = mean(se_p_loo),
                   m_p_waic = mean(p_waic),
                   m_se_waic = mean(se_waic)) 

m_WAIC
#load test results
AllDat_newY <- read.csv(paste0(modelPath, "/rlt/AllDat_newY_",modelversion,".csv"))
AllDat_newY$WMSize <- as.factor(AllDat_newY$WMSize)
levels(AllDat_newY$WMSize) = c("low", "medium",  "high")
AllDat_newY[which(AllDat_newY$Exp == "Exp4"), "Exp"] = "Exp4a"
AllDat_newY[which(AllDat_newY$Exp == "Exp5"), "Exp"] = "Exp4b"

7 Plot Results

7.1 RP

RP_bias_obs  <- ggplot(data = AllExpData%>% filter(Exp != 'Exp5') %>%  
  dplyr::group_by(Exp, curDur, WMSize, NSub) %>% 
  dplyr::summarize(n = n(),
                   m_repDur = mean(repDur),
                   se_repDur = sd(repDur)/sqrt(n-1)) %>%
  dplyr::group_by(Exp, curDur, WMSize) %>% 
  dplyr::summarize(m_m_repDur = mean(m_repDur),
                   m_se_repDur = mean(se_repDur)), aes(x = curDur, y = m_m_repDur-curDur,  color=as.factor(WMSize))) +
  geom_point()+
  geom_line()+
  #geom_errorbar(width=.2,  aes(ymin = m_m_repDur-curDur - m_se_repDur, ymax = m_m_repDur -curDur + m_se_repDur)) +
  geom_hline(yintercept = 0, linetype='dashed')+
  facet_grid(cols = vars(Exp)) +
  labs(x="Sample intervals (s)", y="Reproduction bias(s)", shape=" ", color = "Memory Load")+
  theme_new+colorSet3+guides(shape="none")+
   scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp', 'curDur'. You can override using the `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_obs.png"), RP_bias_obs, width = 6, height = 4)

RP_bias_obs

RP_bias_obs  <- ggplot(data = AllExpData%>% filter(Exp != 'Exp5') %>%  
  dplyr::group_by(Exp, curDur, WMSize, NSub) %>% 
  dplyr::summarize(n = n(),
                   m_repDur = mean(repDur),
                   se_repDur = sd(repDur)/sqrt(n-1)) %>%
  dplyr::group_by(Exp, curDur, WMSize) %>% 
  dplyr::summarize(m_m_repDur = mean(m_repDur),
                   m_se_repDur = mean(se_repDur)), aes(x = curDur, y = m_m_repDur-curDur,  color=as.factor(WMSize))) +
  geom_point()+
  geom_line()+
  geom_errorbar(width=.2,  aes(ymin = m_m_repDur-curDur - m_se_repDur, ymax = m_m_repDur -curDur + m_se_repDur)) +
  geom_hline(yintercept = 0, linetype='dashed')+
  facet_grid(cols = vars(Exp)) +
  labs(x="Sample intervals (s)", y="Reproduction bias(s)", shape=" ", color = "Memory Load")+
  theme_new+colorSet3+guides(shape="none")+
   scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp', 'curDur'. You can override using the `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_obs.png"), RP_bias_obs, width = 6, height = 4)

RP_bias_obs

RP  <- ggplot(data = m_predY, aes(x = curDur, y = m_m_repDur,  color=WMSize, shape = as.factor('Observation'))) +
  geom_point(size=2, alpha = 0.5)+
  geom_line(data= m_newY, aes(x=curDur, y=m_mu_r, color=WMSize)) +
  geom_abline(slope=1, intercept=0)+
  facet_grid(cols = vars(Exp)) +
  labs(x="Sample intervals (s)", y="Reproduction (s)", shape=" ", color = "Memory Load")+
  theme_new+colorSet3+guides(shape="none")+theme_new +theme(legend.position="top")

ggsave(paste0(getwd(), "/", modelPath, "/figures/RP.png"), RP, width = 6, height = 6)

RP

7.2 RP Slope

# calculate the slope of the cv curve
RPslope_model <- function(df) {
  lm(m_repDur ~ log(curDur), data = df)
}

RPslopes <- mpredY_sub %>% 
  dplyr::group_by(NSub, Exp, WMSize) %>% nest() %>%
  mutate(model = map(data, RPslope_model)) %>% 
  mutate(slope = map(model, broom::tidy)) %>%  # get estimates
  unnest(slope, .drop = TRUE) %>% # remove raw data
  select(-std.error,-statistic, -p.value) %>%  # remove unnessary columns
  spread(term, estimate) %>%   # spread stimates
  dplyr::rename(Intercept = `(Intercept)`, slope = `log(curDur)`)  
RPslopes$data <- NULL
RPslopes$model <- NULL
mRPslopes <- RPslopes%>% dplyr::group_by(WMSize, Exp) %>%
  dplyr::summarize(m_Intercept = mean(Intercept), 
                   m_slope = mean(slope),
                   n = n(), 
                   se_slope = sd(slope)/sqrt(n-1),
                   se_Intercept = sd(Intercept)/sqrt(n-1))
## `summarise()` has grouped output by 'WMSize'. You can override using the
## `.groups` argument.
plt_CVslope <- ggplot(mRPslopes, aes(Exp, m_slope, ymin = m_slope - se_slope, ymax = m_slope + se_slope, group =interaction(Exp, WMSize), color = factor(WMSize)), shape = factor(WMSize))+ 
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.2,  position = position_dodge(width = 0.2)) +
  colorSet5+ 
  labs(x = "", y = TeX("Slope of RP"), color = 'Memory Load') +
  theme_new +theme(legend.position="top")

plt_CVslope
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

7.2.1 Anova analysis on slopes of RP

RPslopes$WMSize = as.factor(RPslopes$WMSize)
ezANOVA(data = RPslopes, dv= slope, wid=NSub, within=.(WMSize), between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## $ANOVA
##       Effect DFn DFd         F            p p<.05        ges
## 2        Exp   4  75  1.916369 1.164556e-01       0.08919332
## 3     WMSize   2 150 26.957878 9.954868e-11     * 0.01482364
## 4 Exp:WMSize   8 150  6.842162 1.194810e-07     * 0.01504611
## 
## $`Mauchly's Test for Sphericity`
##       Effect        W          p p<.05
## 3     WMSize 0.937441 0.09160648      
## 4 Exp:WMSize 0.937441 0.09160648      
## 
## $`Sphericity Corrections`
##       Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 3     WMSize 0.9411242 3.125487e-10         * 0.9645461 1.982414e-10         *
## 4 Exp:WMSize 0.9411242 2.602981e-07         * 0.9645461 1.909273e-07         *
ezANOVA(data = RPslopes %>% filter(Exp =='Exp1'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F         p p<.05         ges
## 2 WMSize   2  30 1.628257 0.2131416       0.002398051
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W         p p<.05
## 2 WMSize 0.8320858 0.2761702      
## 
## $`Sphericity Corrections`
##   Effect       GGe     p[GG] p[GG]<.05       HFe    p[HF] p[HF]<.05
## 2 WMSize 0.8562273 0.2172718           0.9557549 0.214484
ezANOVA(data = RPslopes %>% filter(Exp =='Exp2'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F            p p<.05        ges
## 2 WMSize   2  30 21.20007 1.822755e-06     * 0.07646424
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W         p p<.05
## 2 WMSize 0.8461525 0.3105563      
## 
## $`Sphericity Corrections`
##   Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2 WMSize 0.8666656 7.302883e-06         * 0.9698478 2.493639e-06         *
ezANOVA(data = RPslopes %>% filter(Exp =='Exp3'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd         F         p p<.05        ges
## 2 WMSize   2  30 0.4137346 0.6648914       0.00152841
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W         p p<.05
## 2 WMSize 0.9800567 0.8684769      
## 
## $`Sphericity Corrections`
##   Effect       GGe     p[GG] p[GG]<.05      HFe     p[HF] p[HF]<.05
## 2 WMSize 0.9804466 0.6610089           1.126392 0.6648914
ezANOVA(data = RPslopes %>% filter(Exp =='Exp4a'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F           p p<.05        ges
## 2 WMSize   2  30 7.246942 0.002705914     * 0.01861911
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W         p p<.05
## 2 WMSize 0.9730381 0.8258648      
## 
## $`Sphericity Corrections`
##   Effect      GGe       p[GG] p[GG]<.05      HFe       p[HF] p[HF]<.05
## 2 WMSize 0.973746 0.002974256         * 1.117022 0.002705914         *
ezANOVA(data = RPslopes %>% filter(Exp =='Exp4b'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F            p p<.05        ges
## 2 WMSize   2  30 11.66985 0.0001782626     * 0.02763263
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W         p p<.05
## 2 WMSize 0.8407011 0.2968187      
## 
## $`Sphericity Corrections`
##   Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2 WMSize 0.8625903 0.0004111601         * 0.9643404 0.0002213399         *

7.2.2 Anova analysis on Intercept of RP

ezANOVA(data = RPslopes, dv= Intercept, wid=NSub, within=.(WMSize), between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## $ANOVA
##       Effect DFn DFd        F            p p<.05         ges
## 2        Exp   4  75 1.100218 3.628149e-01       0.050380425
## 3     WMSize   2 150 4.729965 1.018536e-02     * 0.006009338
## 4 Exp:WMSize   8 150 6.814933 1.283386e-07     * 0.033669273
## 
## $`Mauchly's Test for Sphericity`
##       Effect         W            p p<.05
## 3     WMSize 0.7988765 0.0002464606     *
## 4 Exp:WMSize 0.7988765 0.0002464606     *
## 
## $`Sphericity Corrections`
##       Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 3     WMSize 0.8325539 1.506819e-02         * 0.8491741 1.449261e-02         *
## 4 Exp:WMSize 0.8325539 1.166472e-06         * 0.8491741 9.364511e-07         *

7.3 RP Bias

  RP_bias  <- ggplot(data = m_predY, aes(x = curDur, y = m_m_repDur - curDur,color=WMSize, shape = as.factor('Observation'))) +
  geom_point(size=2, alpha = 0.5)+
  geom_line(data= m_newY, aes(x=curDur, y=m_mu_r-curDur, color=WMSize)) +
  geom_hline(yintercept = 0, linetype='dashed')+
  facet_grid(cols = vars(Exp)) +
  labs(x=" ", y="Reproduction bias (s)", shape=" ", color = "Memory Load")+theme_new+
  colorSet3+guides(shape="none")

ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias.png"), RP_bias, width = 6, height = 6)

RP_bias

7.4 CV

curDurItem <- unique(m_predY$curDur)
RP_CV <- ggplot(data= m_predY, aes(x=curDur, y= cv, color=WMSize, shape = as.factor('Observation'))) +
  geom_point(size=2, alpha = 0.5)+
  geom_line(data = m_newY, aes(x=curDur, y= m_sig_r/m_mu_r, color=WMSize)) +
  facet_grid(~Exp) +
  labs(x="Duration (s)", y="CV", shape=" ", color = "Memory Load")+ theme_new+
  colorSet3+guides(shape="none")+theme(strip.text.x = element_blank())
RP_CV

ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_CV.png"), RP_CV, width = 6, height = 6)

7.5 CV slope

# calculate the slope of the cv curve
cvSlope_model <- function(df) {
  lm(log(cv_obs) ~ log(curDur), data = df)
}

mpredY <- dplyr::group_by(AllDat_predY, curDur, WMSize, Exp, NSub) %>%
  dplyr::summarize(m_repDur = mean(repDur), 
                   sd_repDur = sd(repDur),
                   n = n(), 
                   sd_repDur = sd(repDur),
                   m_mu_r = mean(mu_r),
                   m_sig_r = mean(sig_r),
                   wp = mean(wp),
                   log_lik =mean(log_lik))
## `summarise()` has grouped output by 'curDur', 'WMSize', 'Exp'. You can override
## using the `.groups` argument.
mpredY$cv_obs <- mpredY$sd_repDur/mpredY$m_repDur

CVslopes <- mpredY %>% 
  dplyr::group_by(NSub, Exp, WMSize) %>% nest()  %>%  # nested data
  mutate(model = map(data, cvSlope_model)) %>%  # linear regression
  mutate(slope = map(model, broom::tidy)) %>%  # get estimates
  unnest(slope, .drop = TRUE) %>% # remove raw data
  select(-std.error,-statistic, -p.value) %>%  # remove unnessary columns
  spread(term, estimate) %>%   # spread stimates
  dplyr::rename(obs_cv_Intercept = `(Intercept)`, obs_cv_slope = `log(curDur)`)  # rename columns
CVslopes$data = NULL
CVslopes$model = NULL
#change the table struction of slopes for spss
CVslopes_list <-split(CVslopes, CVslopes$WMSize) 
CVslopes_spss = NULL
for (i in 1: length(CVslopes_list)){
  temp = CVslopes_list[[i]]
  WMSize = unique(temp$WMSize)
  temp$WMSize = NULL
  colnames(temp) = c('Exp',  'NSub', paste0('obs_cv_Intercept_',WMSize), paste0('obs_cv_slope_',WMSize))
  if(i == 1)
    CVslopes_spss = temp
  else
    CVslopes_spss = left_join(CVslopes_spss, temp, by=c("Exp",  "NSub"))
}
mCVslopes <- CVslopes%>% dplyr::group_by(WMSize, Exp) %>%
  dplyr::summarize(n =n(),
                   m_obs_cv_Intercept = mean(obs_cv_Intercept), 
                   m_obs_cv_slope = mean(obs_cv_slope),
                   se_obs_CV_slope = sd(obs_cv_slope)/sqrt(n-1),
                   se_obs_CV_Intercept = sd(obs_cv_Intercept)/sqrt(n-1))
## `summarise()` has grouped output by 'WMSize'. You can override using the
## `.groups` argument.
plt_CVslope <- ggplot(mCVslopes, aes(Exp, m_obs_cv_slope, ymin = m_obs_cv_slope - se_obs_CV_slope, ymax = m_obs_cv_slope + se_obs_CV_slope, group =interaction(Exp, WMSize), color = factor(WMSize)), shape = factor(WMSize))+ 
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.2,  position = position_dodge(width = 0.2)) +
  colorSet5+ 
  labs(x = "", y = TeX("Slope of CV"), color = 'Memory Load') +
  theme_new +theme(legend.position="top")

plt_CVslope
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_CVslope.png"), plt_CVslope, width = 4, height = 4)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_CVIntecept<- ggplot(mCVslopes, aes(Exp, m_obs_cv_Intercept, ymin = m_obs_cv_Intercept - se_obs_CV_Intercept, ymax = m_obs_cv_Intercept + se_obs_CV_Intercept, group =interaction(Exp, WMSize), color = factor(WMSize)), shape = factor(WMSize))+ 
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.2,  position = position_dodge(width = 0.2)) +
  colorSet5+ 
  labs(x = "", y = TeX("Intercept of CV"), color = 'Memory Load') +
  theme_new +theme(legend.position="top")

plt_CVIntecept
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_Intercept.png"), plt_CVIntecept, width = 4, height = 4)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## Figures in the MS
fig_CV_slope_Intecept<-ggarrange(plt_CVslope, plt_CVIntecept, common.legend = TRUE, ncol=2, nrow=1,  labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_CV_slope_Intecept.png"), fig_CV_slope_Intecept, width = 8, height = 4)
fig_CV_slope_Intecept

## Figures in the MS
fig3<-ggarrange(RP_bias, RP_CV, common.legend = TRUE, ncol=1, nrow=2,  labels = c("a", "b"))
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig3.png"), fig3, width = 6, height = 5)
fig3

8 Indifference Pointtemp_data = sumexpSub_diff %>%

8.1 Indifference Point(linear regression)

#Observed Indifference Point
obs_model <- function(df) {
  lm(repDur ~ curDur, data = df)
}

obs_InP <- AllExpData %>% #filter(curDur %in% c(1.1,1.4)) %>%
  dplyr::group_by(NSub, Exp, WMSize) %>% nest()  %>%
  mutate(model = map(data, obs_model)) %>%  # linear regression
  mutate(slope = map(model, broom::tidy)) %>%  # get estimates
  unnest(slope, .drop = TRUE) %>% # remove raw data
  select(-std.error,-statistic, -p.value) %>%  # remove unnessary columns
  spread(term, estimate) %>%   # spread stimates
  dplyr::rename(Intercept = `(Intercept)`, slope = curDur)  # rename columns
obs_InP$model = NULL
obs_InP$data = NULL
obs_InP$inP = obs_InP$Intercept /(1-obs_InP$slope)
write.csv(obs_InP, file = paste0(modelPath, "/rlt/obs_InP.csv"))
mobs_InP = obs_InP %>% group_by(Exp, WMSize)%>%
  dplyr::summarise(n=n(),
                   m_Intercept = mean(Intercept),
                   se_Intercept= sd(Intercept)/sqrt(n-1),
                   m_inP = mean(inP),
                   se_inP = sd(inP)/sqrt(n-1),
                   m_slope = mean(slope),
                   se_slope = sd(slope)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
plt_InP_linear<- ggplot(data = mobs_InP, aes(x= Exp, y=m_inP, color = WMSize))+
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.3,  aes(ymin = m_inP - se_inP, ymax = m_inP + se_inP), position = position_dodge(width = 0.2)) +theme_new+
  labs(colour = "Memory Load")+colorSet3+
  xlab(' ')+ylab("indifference point (s)")+guides(shape="none")+
  theme(legend.position = "top")


ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_InP_linear.png"), plt_InP_linear, width = 5, height = 5)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_InP_linear
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

RP_Exp2  <- ggplot(data = m_predY%>% filter(Exp == 'Exp2'), aes(x = curDur, y = m_m_repDur,  color=WMSize, shape = as.factor('Observation'))) +
  geom_point(size=2, alpha = 0.5)+
  geom_segment(data =mobs_InP%>%filter(Exp == 'Exp2'), aes(x = 0.3, y = 0.3*m_slope+m_Intercept, xend = 2.2, yend = 2.2*m_slope+m_Intercept), arrow = NULL)+  ##dotted lines
  geom_point(data =mobs_InP%>% filter(Exp == 'Exp2'), aes(x = m_inP, y = 0.1, color=WMSize), shape=21) + ## Intersection 
  geom_segment(data =mobs_InP%>% filter(Exp == 'Exp2'), aes(x=m_inP, y = 0.1, xend = m_inP, yend = m_inP), arrow = NULL)+ ##vertical lines for Intersection
 #geom_line(data= m_newY%>% filter(Exp == 'Exp2'), aes(x=curDur, y=m_mu_r, color=WMSize)) +
  geom_abline(slope=1, intercept=0, linetype = 2)+
  labs(x="Sample intervals (s)", y="Reproduction (s)", shape=" ", color = "Memory Load")+
  theme_new+colorSet3+guides(shape="none")+ ylim(0,2)

ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_Exp2.png"), RP_Exp2, width = 4, height = 4)

RP_Exp2

8.2 Indifference Point (bootstraps)

#### Estimated Indifference Point
mRep_model <- function(df) {
  lm(mu_r ~ curDur, data = df)
}

#Observed Indifference Point
obs_model <- function(df) {
  lm(repDur ~ curDur, data = df)
}

# Custom function to find predicted indifference point 
getPredInP_boot <- function(df, idx){
  vars <- c('NSub', 'Exp', 'WMSize')
  gp_vars = syms(vars) 
  slopes <- df[idx, ] %>% 
    dplyr::group_by(!!!gp_vars) %>% nest()  %>%  # nested data
    mutate(model = map(data, mRep_model)) %>%  # linear regression
    mutate(slope = map(model, broom::tidy)) %>%  # get estimates out
    unnest(slope, .drop = TRUE) %>% # remove raw data
    select(-std.error,-statistic, -p.value) %>%  # remove unnessary clumns
    spread(term, estimate) %>%   # spread stimates
    dplyr::rename(minRP = `(Intercept)`, slope = curDur)  # rename columns
  slopes$inP = slopes$minRP /(1-slopes$slope)
  return(slopes$inP)
}


# Custom function to find observed indifference point 
getRPInP_boot <- function(df, idx){
  vars <- c('NSub', 'Exp', 'WMSize')
  gp_vars = syms(vars) 
  slopes <- df[idx, ] %>% 
    dplyr::group_by(!!!gp_vars) %>% nest()  %>%  # nested data
    mutate(model = map(data, obs_model)) %>%  # linear regression
    mutate(slope = map(model, broom::tidy)) %>%  # get estimates out
    unnest(slope, .drop = TRUE) %>% # remove raw data
    select(-std.error,-statistic, -p.value) %>%  # remove unnessary clumns
    spread(term, estimate) %>%   # spread stimates
    dplyr::rename(minRP = `(Intercept)`, slope = curDur)  # rename columns
  slopes$inP = slopes$minRP /(1-slopes$slope)
  return(slopes$inP)
}
#calculate the bootstrapped 95% confidence intervals 
generateCI = FALSE # tag for generation CI
if(generateCI){
  cilist <- NULL
  for(expname in unique(AllDat_predY$Exp)){
    for(nsub in unique(AllDat_predY$NSub)){
      for(WMsize in unique(AllDat_predY$WMSize)){
        dat = AllDat_predY %>% filter(Exp == expname, NSub == nsub, WMSize ==WMsize)
        set.seed(100) 
        num = 1000
        bs_predInP <- boot(dat, getPredInP_boot, R = num) 
        bs_RPInP <- boot(dat, getRPInP_boot, R = num) 
        ci = data.frame(
          Exp = expname,
          NSub = nsub,
          WMSize = WMsize,
          sd_predInP_boot = sd(bs_predInP$t),
          m_predInP_boot = median(bs_predInP$t),
          sd_RPInP_boot = sd(bs_RPInP$t),
          mRP_InP_boot = median(bs_RPInP$t)
        )
        cilist = data.frame(rbind(cilist, ci))
      }
    }
  }
  write.csv(cilist, file = paste0("ci_list_median_1000.csv"))
}
# load the generated indifference point values and mark the outlier
cilist = read.csv(paste0("ci_list_median_1000.csv"))
cilist$Exp = as.factor(cilist$Exp)
cilist$WMSize= as.factor(cilist$WMSize)
cilist$inPOutlier = FALSE
cilist[which(cilist$mRP_InP_boot > 1.7 |cilist$mRP_InP_boot < 0.5 | cilist$m_predInP_boot < 0.5|cilist$m_predInP_boot > 1.7),"inPOutlier"] = TRUE

#check if the outlier is the same as the outliers in variable slope_pr
cilist%>% filter(inPOutlier == TRUE)

8.3 Inifference Point (curve)

AllDat_newY$predErr = AllDat_newY$mu_r - AllDat_newY$curDur

temp_newY <- AllDat_newY %>% filter(curDur > 0.8, curDur < 1.1) %>% select(Exp, WMSize, NSub, predErr, curDur)

InP_curve<- temp_newY%>% dplyr::group_by(Exp, WMSize, NSub)%>%
  dplyr::summarise(minErr = min(abs(predErr)), idx = which.min(abs(predErr)))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
InP_curve$InP_curve = temp_newY[InP_curve$idx,]$curDur 
InP_curve$y = temp_newY[InP_curve$idx,]$predErr + temp_newY[InP_curve$idx,]$curDur
InP_curve

8.4 Inifference Point (Equation)

#sig_sm2 =  sig_s2 + ls*(size[i]);
#wp_pr[i] = sig_sm2 /(sig_sm2 + sigma_pr2); 

Bayparlist$wp_1 = (Bayparlist$sig_s2 + Bayparlist$ls *1 )/ (Bayparlist$sig_s2 + Bayparlist$ls *1 + Bayparlist$sig_pr2_log)

Bayparlist$wp_3 = (Bayparlist$sig_s2 + Bayparlist$ls *2) / (Bayparlist$sig_s2 + Bayparlist$ls *2 + Bayparlist$sig_pr2_log)

Bayparlist$wp_5 = (Bayparlist$sig_s2 + Bayparlist$ls *3) / (Bayparlist$sig_s2 + Bayparlist$ls *3 + Bayparlist$sig_pr2_log)

#log(IP)=[kr*M + (1-w_p)ks*M+ w_p*mu_p]/w_p 

size = 1 
Bayparlist$InP_1 = exp((Bayparlist$kr*size + (1-Bayparlist$wp_1)*Bayparlist$ks*size + Bayparlist$wp_1* Bayparlist$mu_pr_log)/Bayparlist$wp_1)
size = 2
Bayparlist$InP_3 = exp((Bayparlist$kr*size + (1-Bayparlist$wp_3)*Bayparlist$ks*size + Bayparlist$wp_3* Bayparlist$mu_pr_log)/Bayparlist$wp_3)
size = 3
Bayparlist$InP_5 = exp((Bayparlist$kr*size + (1-Bayparlist$wp_5)*Bayparlist$ks*size + Bayparlist$wp_5* Bayparlist$mu_pr_log)/Bayparlist$wp_5)
Bayparlist$sig2_post_1 = (Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*1))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*1) *0.5

Bayparlist$sig2_post_3 =(Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*3))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*3) *0.5 

Bayparlist$sig2_post_5 = (Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*5))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*5) *0.5

Bayparlist%>%filter(Exp == 'Exp2')%>%select("NSub", "Exp", "sig2_post_1", "sig2_post_3","sig2_post_5")
Bayparlist$InP_old_1 = exp((Bayparlist$kr*size - (1-Bayparlist$wp_1)*Bayparlist$ks*size +(Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*1))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*1) *0.5 + Bayparlist$wp_1* Bayparlist$mu_pr_log)/Bayparlist$wp_1)

Bayparlist$InP_old_3 = exp((Bayparlist$kr*size - (1-Bayparlist$wp_3)*Bayparlist$ks*size +(Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*3))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*3) *0.5 + Bayparlist$wp_3* Bayparlist$mu_pr_log)/Bayparlist$wp_3)

Bayparlist$InP_old_5 = exp((Bayparlist$kr*size - (1-Bayparlist$wp_5)*Bayparlist$ks*size +(Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*5))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*5) *0.5 + Bayparlist$wp_5* Bayparlist$mu_pr_log)/Bayparlist$wp_5)
Bayparlist%>%select(NSub, Exp, wp_1, wp_3, wp_5, InP_1, InP_3, InP_5, InP_old_1, InP_old_3, InP_old_5)%>%
  unite(newcol1,  wp_1, InP_1, InP_old_1)%>%
  unite(newcol3,  wp_3, InP_3, InP_old_3)%>%
  unite(newcol5, wp_5, InP_5, InP_old_5)%>%
  melt(id.vars = c("NSub", "Exp")) %>%
  dplyr::rename(
    WMSize = variable,
    newcol = value
  )%>%
  separate(newcol,  c("wp", "InP_Eq", "InP_Eq_old"), sep = "_") -> Baypar_WMSize
Baypar_WMSize$WMSize <- substr(Baypar_WMSize$WMSize, 7, 7)
Baypar_WMSize$WMSize <- as.factor(Baypar_WMSize$WMSize)
Baypar_WMSize$InP_Eq = as.numeric(Baypar_WMSize$InP_Eq)
Baypar_WMSize$InP_Eq_old = as.numeric(Baypar_WMSize$InP_Eq_old)
levels(Baypar_WMSize$WMSize) = c("low", "medium",  "high")


##save csv for SPSS
left_join(Baypar_WMSize, InP_curve%>%select("Exp", "WMSize","NSub", "InP_curve"),  by = c("Exp","WMSize","NSub")) -> Baypar_WMSize

left_join(Baypar_WMSize, cilist %>% select("Exp", "NSub", "WMSize", "m_predInP_boot", "mRP_InP_boot"), by = c("Exp","WMSize","NSub"))-> Baypar_WMSize

8.5 check indifference Points

Baypar_WMSize_melt <- melt(Baypar_WMSize%>%select("Exp","WMSize","NSub", "InP_Eq", "InP_Eq_old", "InP_curve", "m_predInP_boot", "mRP_InP_boot"), id.vars = c("Exp","WMSize","NSub"),
                           variable.name = "Type",
                           value.name = "InP" )


Baypar_WMSize_melt %>%filter(Type== "InP_curve") %>% dplyr::group_by(Exp)%>%
  dplyr::summarise(n = n(),
                   m_InP = mean(InP), se_InP= sd(InP)/sqrt(n-1))
Baypar_WMSize_melt$Exp = as.factor(Baypar_WMSize_melt$Exp)

Baypar_WMSize_melt%>% dplyr::group_by(Exp, WMSize, Type)%>%
  dplyr::summarise(n = n(),
                   m_InP = mean(InP), se_InP= sd(InP)/sqrt(n-1))-> mBaypar_InP
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
#plot InP_curve(the intersections of the Prediction curve with the diagonal)
plt_InP<- ggplot(data = mBaypar_InP%>%filter(Type== "InP_curve"), aes(x= Exp, y=m_InP, color = WMSize, shape = Type))+
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.2,  aes(ymin = m_InP - se_InP, ymax = m_InP + se_InP), position = position_dodge(width = 0.2)) +theme_new+
  labs(colour = "Memory Load")+colorSet3+
  xlab(' ')+ylab("indifference point (s)")+guides(shape="none")+
  theme(legend.position = "top")


ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_InP.png"), plt_InP, width = 3, height = 3)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_InP
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

temp_inp = mBaypar_InP%>%filter(Type== "InP_curve", Exp == 'Exp2')
RP_bias_Exp2  <- ggplot(data = m_predY%>%filter(Exp == 'Exp2'), aes(x = curDur, y = m_m_repDur - curDur,color=WMSize, shape = as.factor('Observation'))) +
  geom_point(size=2, alpha = 0.5)+
  geom_line(data= m_newY%>%filter(Exp == 'Exp2'), aes(x=curDur, y=m_mu_r-curDur, color=WMSize)) +
  geom_hline(yintercept = 0, linetype='dashed')+
  geom_point(data =temp_inp, aes(x = m_InP+0.03, y = -0.5, color=WMSize), shape=21) + ## Intersection 
  geom_segment(data =temp_inp, aes(x=m_InP+0.03, y = -0.5, xend = m_InP+0.03, yend = 0), arrow = NULL)+ ##vertical lines for Intersection
  labs(x=" ", y="Reproduction bias (s)", shape=" ", color = "Memory Load")+theme_new+
  colorSet3+guides(shape="none")

ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_Exp2.png"), RP_bias_Exp2, width = 5, height = 3)

RP_bias_Exp2

8.6 anova on observed InP

ezANOVA(data = Baypar_WMSize, dv= InP_curve, wid=NSub, within = .(WMSize), between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## $ANOVA
##       Effect DFn DFd        F            p p<.05         ges
## 2        Exp   4  75 1.198179 3.187005e-01       0.054696877
## 3     WMSize   2 150 2.377013 9.631587e-02       0.002987228
## 4 Exp:WMSize   8 150 4.618376 4.655267e-05     * 0.022755616
## 
## $`Mauchly's Test for Sphericity`
##       Effect         W            p p<.05
## 3     WMSize 0.3471828 1.001507e-17     *
## 4 Exp:WMSize 0.3471828 1.001507e-17     *
## 
## $`Sphericity Corrections`
##       Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 3     WMSize 0.6050276 0.1210006309           0.6095966 0.1207130064          
## 4 Exp:WMSize 0.6050276 0.0009604633         * 0.6095966 0.0009269622         *
pairwise.t.test(Baypar_WMSize$InP_curve, Baypar_WMSize$Exp)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  Baypar_WMSize$InP_curve and Baypar_WMSize$Exp 
## 
##       Exp1  Exp2  Exp3  Exp4a
## Exp2  1.000 -     -     -    
## Exp3  1.000 1.000 -     -    
## Exp4a 1.000 0.951 1.000 -    
## Exp4b 0.075 0.019 0.031 0.604
## 
## P value adjustment method: holm
ezANOVA(data = Baypar_WMSize %>% filter(Exp =='Exp2'), dv= InP_curve, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F           p p<.05       ges
## 2 WMSize   2  30 5.834043 0.007240452     * 0.0340583
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W            p p<.05
## 2 WMSize 0.1674321 3.688675e-06     *
## 
## $`Sphericity Corrections`
##   Effect       GGe      p[GG] p[GG]<.05       HFe      p[HF] p[HF]<.05
## 2 WMSize 0.5456824 0.02548011         * 0.5558358 0.02476867         *
ezANOVA(data = Baypar_WMSize %>% filter(Exp =='Exp3'), dv= InP_curve, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F           p p<.05        ges
## 2 WMSize   2  30 8.753035 0.001012937     * 0.03895252
## 
## $`Mauchly's Test for Sphericity`
##   Effect          W            p p<.05
## 2 WMSize 0.06186308 3.467534e-09     *
## 
## $`Sphericity Corrections`
##   Effect       GGe       p[GG] p[GG]<.05       HFe       p[HF] p[HF]<.05
## 2 WMSize 0.5159594 0.009073344         * 0.5194236 0.008930414         *
ezANOVA(data = Baypar_WMSize %>% filter(Exp =='Exp4a'), dv= InP_curve, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd         F         p p<.05        ges
## 2 WMSize   2  30 0.7386935 0.4862266       0.00420096
## 
## $`Mauchly's Test for Sphericity`
##   Effect          W            p p<.05
## 2 WMSize 0.07926534 1.965999e-08     *
## 
## $`Sphericity Corrections`
##   Effect       GGe    p[GG] p[GG]<.05       HFe     p[HF] p[HF]<.05
## 2 WMSize 0.5206341 0.408452           0.5251298 0.4094818
ezANOVA(data = Baypar_WMSize %>% filter(Exp =='Exp4b'), dv= InP_curve, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd       F          p p<.05        ges
## 2 WMSize   2  30 3.57389 0.04053459     * 0.04484435
## 
## $`Mauchly's Test for Sphericity`
##   Effect        W           p p<.05
## 2 WMSize 0.450539 0.003768136     *
## 
## $`Sphericity Corrections`
##   Effect       GGe     p[GG] p[GG]<.05       HFe      p[HF] p[HF]<.05
## 2 WMSize 0.6453857 0.0647428           0.6802842 0.06184001

9 Export data for spss

### change the table struction for spss
Bayparlist_list <-split(Baypar_WMSize %>% select(c("NSub","Exp","WMSize", "wp","InP_curve")), Baypar_WMSize$WMSize) 
Bayparlist_spss = NULL
for (i in 1: length(Bayparlist_list)){
  temp = Bayparlist_list[[i]]
  WMSize = unique(temp$WMSize)
  temp$WMSize = NULL
  colnames(temp) = c('NSub','Exp', paste0('wp_',WMSize), paste0('InP_curve_',WMSize))
  if(i == 1)
    Bayparlist_spss = temp
  else
    Bayparlist_spss = left_join(Bayparlist_spss, temp, by=c("Exp",  "NSub"))
}
write_csv(Bayparlist_spss, paste0(modelPath, '/rlt/Bayparlist_spss.csv'))

10 anova analysis

10.1 Anova on mean reproduction biases

mpredY_sub$RP_Bias = mpredY_sub$m_repDur-mpredY_sub$curDur
RP_bias_Anova <- ezANOVA(data = mpredY_sub, dv= RP_Bias, wid=NSub, within= .(curDur, WMSize), between = .(Exp), detailed = TRUE, return_aov = FALSE )
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: "curDur" will be treated as numeric.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Warning: There is at least one numeric within variable, therefore aov() will be
## used for computation and no assumption checks will be obtained.
RP_bias_Anova
## $ANOVA
##              Effect DFn DFd         SSn       SSd          F            p p<.05
## 1               Exp   4  75  0.37301817 6.5861689   1.061936 3.813811e-01      
## 2            curDur   1  75 49.69246885 7.6805406 485.243856 1.753714e-34     *
## 4            WMSize   2 150  0.04636729 0.6908236   5.033914 7.656681e-03     *
## 3        Exp:curDur   4  75  0.78848535 7.6805406   1.924878 1.150304e-01      
## 5        Exp:WMSize   8 150  0.25333639 0.6908236   6.875933 1.093463e-07     *
## 6     curDur:WMSize   2 150  0.12351691 0.3507396  26.412100 1.488886e-10     *
## 7 Exp:curDur:WMSize   8 150  0.12715860 0.3507396   6.797704 1.342826e-07     *
##           ges
## 1 0.023787466
## 2 0.764490798
## 4 0.003019758
## 3 0.048984109
## 5 0.016279575
## 6 0.008004056
## 7 0.008238098
# main effect of Duration  F(1.177, 3.532) = 377.965, p < .001, ηp² = .863.
(RP_bias_Anova$ANOVA)$DFn[3] *(RP_bias_Anova$`Sphericity Corrections`)$GGe[1]
## numeric(0)
(RP_bias_Anova$ANOVA)$DFd[3] *(RP_bias_Anova$`Sphericity Corrections`)$GGe[1]
## numeric(0)
#Duration × Experiment,  F(12, 240) = 2.506, p = .004, ηp² = .111
(RP_bias_Anova$ANOVA)$DFn[5] *(RP_bias_Anova$`Sphericity Corrections`)$GGe[2]
## numeric(0)
(RP_bias_Anova$ANOVA)$DFd[5] *(RP_bias_Anova$`Sphericity Corrections`)$GGe[2]
## numeric(0)
mpredY_sub <- ungroup(mpredY_sub)
res.aov <-  rstatix::anova_test(data = mpredY_sub, dv = RP_Bias, wid = NSub, within = c(curDur, WMSize), between = Exp)
## Warning: The 'wid' column contains duplicate ids across between-subjects
## variables. Automatic unique id will be created
get_anova_table(res.aov, correction = "GG")

10.2 variance of prior (anova)

ezANOVA(data = Bayparlist, dv= sig_pr2_log, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
##   Effect DFn DFd         F         p p<.05        ges
## 1    Exp   4  75 0.7160159 0.5836025       0.03678287
## 
## $`Levene's Test for Homogeneity of Variance`
##   DFn DFd       SSn      SSd         F         p p<.05
## 1   4  75 0.1582023 3.779256 0.7848881 0.5385686

10.3 variance of motor noise (anova)

ezANOVA(data = Bayparlist, dv= sig_mn2, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
##   Effect DFn DFd        F         p p<.05        ges
## 1    Exp   4  75 1.864222 0.1255674       0.09043378
## 
## $`Levene's Test for Homogeneity of Variance`
##   DFn DFd          SSn         SSd        F         p p<.05
## 1   4  75 0.0005070715 0.007658422 1.241456 0.3007183

10.4 Ks anova

ezANOVA(data = Bayparlist, dv= ks, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
##   Effect DFn DFd        F            p p<.05       ges
## 1    Exp   4  75 21.71153 6.283521e-12     * 0.5365969
## 
## $`Levene's Test for Homogeneity of Variance`
##   DFn DFd       SSn      SSd        F            p p<.05
## 1   4  75 0.4637314 1.347294 6.453651 0.0001614605     *

10.5 mean of prior (anova)

ezANOVA(data = Bayparlist, dv= mu_pr, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
##   Effect DFn DFd        F         p p<.05        ges
## 1    Exp   4  75 1.372231 0.2516599       0.06819476
## 
## $`Levene's Test for Homogeneity of Variance`
##   DFn DFd       SSn      SSd         F         p p<.05
## 1   4  75 0.8682032 27.17694 0.5989934 0.6644974

11 weight of prior

11.1 Weight of the prior \(w_p\)

Baypar_WMSize$wp = as.numeric(Baypar_WMSize$wp)
mwp <- Baypar_WMSize%>%dplyr::group_by(Exp, WMSize)%>% dplyr::summarise(m_wp = mean(wp), n= n(), se_wp = sd(wp)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
plt_wp <- ggplot(mwp, aes(Exp, m_wp, ymin = m_wp - se_wp, ymax = m_wp + se_wp, group =interaction(Exp, WMSize), color = factor(WMSize)), shape = factor(WMSize))+ 
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.2,  position = position_dodge(width = 0.2)) +
  #coord_cartesian(ylim = c(0.5, 1)) +
  colorSet5+
  labs(x = "", y = TeX("Weight of the prior $w_p$"), color = 'Memory Load') +
  theme_new + theme(legend.position="top")

plt_wp
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_wp.png"), plt_wp, width = 3, height = 3)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

11.2 combine InP and wp

fig4<-ggarrange(plt_wp, plt_InP, common.legend = TRUE, ncol=2, nrow=1,  labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig4.png"), fig4, width = 6, height = 3)
fig4

fig5<-ggarrange(RP_bias_obs, plt_InP_linear, common.legend = TRUE, ncol=2, nrow=1, widths = c(5,3), labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig5.png"), fig5, width = 6, height = 3)
fig5

11.3 ANOVA analysis on wp

ezANOVA(data = Baypar_WMSize, dv= wp, wid=NSub, within=.(WMSize), between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## $ANOVA
##       Effect DFn DFd          F            p p<.05         ges
## 2        Exp   4  75   2.249548 7.163591e-02       0.106466774
## 3     WMSize   2 150 123.376009 2.079259e-32     * 0.011162610
## 4 Exp:WMSize   8 150  25.074955 3.255859e-24     * 0.009093747
## 
## $`Mauchly's Test for Sphericity`
##       Effect         W            p p<.05
## 3     WMSize 0.0228096 1.779717e-61     *
## 4 Exp:WMSize 0.0228096 1.779717e-61     *
## 
## $`Sphericity Corrections`
##       Effect       GGe        p[GG] p[GG]<.05      HFe        p[HF] p[HF]<.05
## 3     WMSize 0.5057682 1.120300e-17         * 0.506003 1.102351e-17         *
## 4 Exp:WMSize 0.5057682 2.492781e-13         * 0.506003 2.463140e-13         *
ezANOVA(data = Baypar_WMSize%>%filter(Exp =='Exp2'), dv= wp, wid=NSub, within= .(WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd       F            p p<.05        ges
## 2 WMSize   2  30 55.0847 9.057763e-11     * 0.04339272
## 
## $`Mauchly's Test for Sphericity`
##   Effect          W            p p<.05
## 2 WMSize 0.02641187 8.965915e-12     *
## 
## $`Sphericity Corrections`
##   Effect       GGe        p[GG] p[GG]<.05      HFe        p[HF] p[HF]<.05
## 2 WMSize 0.5066913 1.869688e-06         * 0.508133 1.815792e-06         *
ezANOVA(data = Baypar_WMSize%>%filter(Exp =='Exp4a'), dv= wp, wid=NSub, within= .(WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F            p p<.05        ges
## 2 WMSize   2  30 33.98248 1.953251e-08     * 0.01772593
## 
## $`Mauchly's Test for Sphericity`
##   Effect          W           p p<.05
## 2 WMSize 0.01054139 1.44639e-14     *
## 
## $`Sphericity Corrections`
##   Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2 WMSize 0.5026493 3.187613e-05         * 0.5032182 3.160504e-05         *
ezANOVA(data = Baypar_WMSize%>%filter(Exp =='Exp4b'), dv= wp, wid=NSub, within= .(WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##   Effect DFn DFd        F            p p<.05        ges
## 2 WMSize   2  30 34.92081 1.469415e-08     * 0.02125106
## 
## $`Mauchly's Test for Sphericity`
##   Effect          W            p p<.05
## 2 WMSize 0.02260196 3.013161e-12     *
## 
## $`Sphericity Corrections`
##   Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2 WMSize 0.5057151 2.627002e-05         * 0.5069454 2.578052e-05         *
### pairwise.t.test on wp
pairwise.t.test(Baypar_WMSize$wp, Baypar_WMSize$Exp)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  Baypar_WMSize$wp and Baypar_WMSize$Exp 
## 
##       Exp1    Exp2    Exp3    Exp4a  
## Exp2  0.00016 -       -       -      
## Exp3  0.12532 0.21889 -       -      
## Exp4a 1.00000 0.00029 0.16667 -      
## Exp4b 1.00000 0.00323 0.47867 1.00000
## 
## P value adjustment method: holm
##independent T test
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp1')))$wp)
## 
##  One Sample t-test
## 
## data:  (Baypar_WMSize %>% filter(Exp %in% c("Exp1")))$wp
## t = 14.032, df = 47, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.3459136 0.4617008
## sample estimates:
## mean of x 
## 0.4038072
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp2')))$wp)
## 
##  One Sample t-test
## 
## data:  (Baypar_WMSize %>% filter(Exp %in% c("Exp2")))$wp
## t = 16.632, df = 47, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.5184761 0.6611631
## sample estimates:
## mean of x 
## 0.5898196
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp3')))$wp)
## 
##  One Sample t-test
## 
## data:  (Baypar_WMSize %>% filter(Exp %in% c("Exp3")))$wp
## t = 20.765, df = 47, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.4554826 0.5532064
## sample estimates:
## mean of x 
## 0.5043445
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp4a')))$wp)
## 
##  One Sample t-test
## 
## data:  (Baypar_WMSize %>% filter(Exp %in% c("Exp4a")))$wp
## t = 14.029, df = 47, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.3520475 0.4699193
## sample estimates:
## mean of x 
## 0.4109834
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp4b')))$wp)
## 
##  One Sample t-test
## 
## data:  (Baypar_WMSize %>% filter(Exp %in% c("Exp4b")))$wp
## t = 14.529, df = 47, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.3777757 0.4992013
## sample estimates:
## mean of x 
## 0.4384885

11.4 standard variance of Ds \(\sigma_{s}\)

ezANOVA(data = Bayparlist, dv= sig_s2, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
##   Effect DFn DFd         F         p p<.05        ges
## 1    Exp   4  75 0.3845193 0.8190605       0.02009558
## 
## $`Levene's Test for Homogeneity of Variance`
##   DFn DFd        SSn       SSd         F         p p<.05
## 1   4  75 0.02897922 0.8327796 0.6524659 0.6269356
### pairwise.t.test on standard variance of $D_s$
pairwise.t.test(Bayparlist$sig_s2, Bayparlist$Exp)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  Bayparlist$sig_s2 and Bayparlist$Exp 
## 
##       Exp1 Exp2 Exp3 Exp4a
## Exp2  1    -    -    -    
## Exp3  1    1    -    -    
## Exp4a 1    1    1    -    
## Exp4b 1    1    1    1    
## 
## P value adjustment method: holm
##independent T test
t.test((Bayparlist%>%filter(Exp %in%c('Exp1')))$sig_s2)
## 
##  One Sample t-test
## 
## data:  (Bayparlist %>% filter(Exp %in% c("Exp1")))$sig_s2
## t = 7.2797, df = 15, p-value = 2.698e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.04681825 0.08558531
## sample estimates:
##  mean of x 
## 0.06620178
t.test((Bayparlist%>%filter(Exp %in%c('Exp2')))$sig_s2)
## 
##  One Sample t-test
## 
## data:  (Bayparlist %>% filter(Exp %in% c("Exp2")))$sig_s2
## t = 2.5569, df = 15, p-value = 0.0219
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.01307008 0.14402351
## sample estimates:
##  mean of x 
## 0.07854679
t.test((Bayparlist%>%filter(Exp %in%c('Exp3')))$sig_s2)
## 
##  One Sample t-test
## 
## data:  (Bayparlist %>% filter(Exp %in% c("Exp3")))$sig_s2
## t = 11.647, df = 15, p-value = 6.497e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.06165666 0.08927951
## sample estimates:
##  mean of x 
## 0.07546809
t.test((Bayparlist%>%filter(Exp %in%c('Exp4a')))$sig_s2)
## 
##  One Sample t-test
## 
## data:  (Bayparlist %>% filter(Exp %in% c("Exp4a")))$sig_s2
## t = 7.0467, df = 15, p-value = 3.959e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.03222067 0.06016474
## sample estimates:
##  mean of x 
## 0.04619271
t.test((Bayparlist%>%filter(Exp %in%c('Exp4b')))$sig_s2)
## 
##  One Sample t-test
## 
## data:  (Bayparlist %>% filter(Exp %in% c("Exp4b")))$sig_s2
## t = 1.7736, df = 15, p-value = 0.09643
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.01870072  0.20404884
## sample estimates:
##  mean of x 
## 0.09267406
# calculate conhens d
cohensD(sig_s2 ~ Exp,
        data   = Bayparlist%>%filter(Exp %in% c('Exp1', 'Exp2')),
        method = "paired")
## Warning in cohensD(sig_s2 ~ Exp, data = Bayparlist %>% filter(Exp %in%
## c("Exp1", : calculating paired samples Cohen's d using formula input. Results
## will be incorrect if cases do not appear in the same order for both levels of
## the grouping factor
## [1] 0.08918868
cohensD(sig_s2 ~ Exp,
        data   = Bayparlist%>%filter(Exp %in% c('Exp2', 'Exp3')),
        method = "paired")
## Warning in cohensD(sig_s2 ~ Exp, data = Bayparlist %>% filter(Exp %in%
## c("Exp2", : calculating paired samples Cohen's d using formula input. Results
## will be incorrect if cases do not appear in the same order for both levels of
## the grouping factor
## [1] 0.02474716
cohensD(sig_s2 ~ Exp,
        data   = Bayparlist%>%filter(Exp %in% c('Exp3', 'Exp4a')),
        method = "paired")
## Warning in cohensD(sig_s2 ~ Exp, data = Bayparlist %>% filter(Exp %in%
## c("Exp3", : calculating paired samples Cohen's d using formula input. Results
## will be incorrect if cases do not appear in the same order for both levels of
## the grouping factor
## [1] 0.8774215
cohensD(sig_s2 ~ Exp,
        data   = Bayparlist%>%filter(Exp %in% c('Exp2', 'Exp4a')),
        method = "paired")
## Warning in cohensD(sig_s2 ~ Exp, data = Bayparlist %>% filter(Exp %in%
## c("Exp2", : calculating paired samples Cohen's d using formula input. Results
## will be incorrect if cases do not appear in the same order for both levels of
## the grouping factor
## [1] 0.2541347

11.5 cv slope

ezANOVA(data = CVslopes, dv= obs_cv_slope, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Warning: Collapsing data to cell means. *IF* the requested effects are a subset
## of the full design, you must use the "within_full" argument, else results may be
## inaccurate.
## Coefficient covariances computed by hccm()
## $ANOVA
##   Effect DFn DFd          F         p p<.05         ges
## 1    Exp   4  75 0.08928208 0.9855437       0.004739144
## 
## $`Levene's Test for Homogeneity of Variance`
##   DFn DFd        SSn      SSd        F         p p<.05
## 1   4  75 0.08788252 1.346571 1.223699 0.3079875
### pairwise.t.test on standard variance of $D_s$
pairwise.t.test(CVslopes$obs_cv_slope, CVslopes$Exp)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  CVslopes$obs_cv_slope and CVslopes$Exp 
## 
##       Exp1 Exp2 Exp3 Exp4a
## Exp2  1    -    -    -    
## Exp3  1    1    -    -    
## Exp4a 1    1    1    -    
## Exp4b 1    1    1    1    
## 
## P value adjustment method: holm
##independent T test
t.test((CVslopes%>%filter(Exp %in%c('Exp1')))$obs_cv_slope)
## 
##  One Sample t-test
## 
## data:  (CVslopes %>% filter(Exp %in% c("Exp1")))$obs_cv_slope
## t = -6.5917, df = 47, p-value = 3.406e-08
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.2970793 -0.1581479
## sample estimates:
##  mean of x 
## -0.2276136
t.test((CVslopes%>%filter(Exp %in%c('Exp2')))$obs_cv_slope)
## 
##  One Sample t-test
## 
## data:  (CVslopes %>% filter(Exp %in% c("Exp2")))$obs_cv_slope
## t = -7.47, df = 47, p-value = 1.59e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.3369816 -0.1939864
## sample estimates:
## mean of x 
## -0.265484
t.test((CVslopes%>%filter(Exp %in%c('Exp3')))$obs_cv_slope)
## 
##  One Sample t-test
## 
## data:  (CVslopes %>% filter(Exp %in% c("Exp3")))$obs_cv_slope
## t = -6.932, df = 47, p-value = 1.037e-08
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.3153826 -0.1735021
## sample estimates:
##  mean of x 
## -0.2444424
t.test((CVslopes%>%filter(Exp %in%c('Exp4a')))$obs_cv_slope)
## 
##  One Sample t-test
## 
## data:  (CVslopes %>% filter(Exp %in% c("Exp4a")))$obs_cv_slope
## t = -5.4143, df = 47, p-value = 2.047e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.3547829 -0.1625578
## sample estimates:
##  mean of x 
## -0.2586703
t.test((CVslopes%>%filter(Exp %in%c('Exp4b')))$obs_cv_slope)
## 
##  One Sample t-test
## 
## data:  (CVslopes %>% filter(Exp %in% c("Exp4b")))$obs_cv_slope
## t = -7.1687, df = 47, p-value = 4.537e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.2997242 -0.1683650
## sample estimates:
##  mean of x 
## -0.2340446
# calculate conhens d
cohensD(obs_cv_slope ~ Exp,
        data   = CVslopes%>%filter(Exp %in% c('Exp1', 'Exp2')),
        method = "paired")
## Warning in cohensD(obs_cv_slope ~ Exp, data = CVslopes %>% filter(Exp %in% :
## calculating paired samples Cohen's d using formula input. Results will be
## incorrect if cases do not appear in the same order for both levels of the
## grouping factor
## [1] 0.1178123
cohensD(obs_cv_slope ~ Exp,
        data   = CVslopes%>%filter(Exp %in% c('Exp2', 'Exp3')),
        method = "paired")
## Warning in cohensD(obs_cv_slope ~ Exp, data = CVslopes %>% filter(Exp %in% :
## calculating paired samples Cohen's d using formula input. Results will be
## incorrect if cases do not appear in the same order for both levels of the
## grouping factor
## [1] 0.06429713
cohensD(obs_cv_slope ~ Exp,
        data   = CVslopes%>%filter(Exp %in% c('Exp3', 'Exp4a')),
        method = "paired")
## Warning in cohensD(obs_cv_slope ~ Exp, data = CVslopes %>% filter(Exp %in% :
## calculating paired samples Cohen's d using formula input. Results will be
## incorrect if cases do not appear in the same order for both levels of the
## grouping factor
## [1] 0.04045837
cohensD(obs_cv_slope ~ Exp,
        data   = CVslopes%>%filter(Exp %in% c('Exp2', 'Exp4a')),
        method = "paired")
## Warning in cohensD(obs_cv_slope ~ Exp, data = CVslopes %>% filter(Exp %in% :
## calculating paired samples Cohen's d using formula input. Results will be
## incorrect if cases do not appear in the same order for both levels of the
## grouping factor
## [1] 0.01404826

12 Model prediction error

m_predErr_sub<- mpredY_sub%>% 
  dplyr::group_by(Exp, WMSize, NSub) %>% dplyr::summarise(
    mpredRP_err=mean(predRP_err),
    mpredVar_err=mean(predVar_err),
    mpredcv_err = mean(predcv_err),
    mpredRP_rerr = mean(predRP_rerr),
    mpredVar_rerr = mean(predVar_rerr),
    mpredcv_rerr = mean(predcv_rerr))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
m_predErr<- m_predY%>% 
  dplyr::group_by(Exp, WMSize) %>% dplyr::summarise(
    mmpredcv_err = mean(mpredcv_err),
    mmpredRP_err=mean(mpredRP_err),
    mmpredVar_err=mean(mpredVar_err),
    mmpredRP_rerr = mean(mpredRP_rerr),
    mmpredVar_rerr = mean(mpredVar_rerr),
    mmpredcv_rerr = mean(mpredcv_rerr))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
m_predErr
plt_rErrorScatter = ggplot(m_predErr_sub, aes(mpredRP_rerr*100, mpredcv_rerr*100, color = WMSize, alpha = .9)) + 
  #geom_hline(yintercept = 0, linetype='dashed')+ geom_vline(xintercept = 0, linetype='dashed')+ 
  geom_point() +
  geom_point(data = m_predErr, aes(mmpredRP_rerr*100, mmpredcv_rerr*100, color = WMSize, alpha = .9, size = 1 ))+
  xlab('Relative prediction error of reproduction mean(%)')+ ylab('Relative prediction error of CV (%)')+colorSet3+
  facet_wrap(~Exp)+
  theme_new+ theme(legend.position = 'top')+guides(size="none")+guides(alpha="none")

plt_rErrorScatter

plt_rErrorScatter = ggplot(m_predErr_sub, aes(mpredRP_rerr*100, mpredVar_rerr*100, color = WMSize, alpha = .9)) + 
  #geom_hline(yintercept = 0, linetype='dashed')+ geom_vline(xintercept = 0, linetype='dashed')+ 
  geom_point() +
  geom_point(data = m_predErr, aes(mmpredRP_rerr*100, mmpredVar_rerr*100, color = WMSize, alpha = .9, size = 1 ))+
  xlab('Relative prediction error in the RP means (%)')+ ylab('Relative prediction error in the RP variance (%)')+colorSet3+
  facet_wrap(~Exp)+
  theme_new+ theme(legend.position = 'top')+guides(size="none")+guides(alpha="none")

plt_rErrorScatter

plt_ErrorScatter = ggplot(m_predErr_sub, aes(mpredRP_err, mpredVar_err, color = WMSize, alpha = .9)) + 
  geom_hline(yintercept = 0, linetype='dashed')+ geom_vline(xintercept = 0, linetype='dashed')+ 
  geom_point() +
  geom_point(data = m_predErr, aes(mmpredRP_err, mmpredVar_err, color = WMSize, alpha = .9, size = 1 ))+
  xlab('Prediction error in the RP means (s)')+ ylab('Prediction error in the RP variance (s)')+colorSet3+
  facet_wrap(~Exp)+
  theme_new+ theme(legend.position = 'top')+guides(size="none")+guides(alpha="none")

ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_ErrorScatter.png"), plt_ErrorScatter, width = 7, height = 4)

plt_ErrorScatter

13 model comparison (logarithmic vs. linear)

m_predErr_sub$model = 'logarithmic'
m_predErr$model = 'logarithmic'
linear_model = 'linear_rstan'
m_predErr_linear = read.csv(paste0(getwd(), "/", rstanmodelPath, '/models/', linear_model, "/rlt/m_predErr_", linear_model, ".csv"))
m_predErr_linear$X = NULL
m_predErr_sub_linear = read.csv(paste0(getwd(), "/", rstanmodelPath, '/models/', linear_model, "/rlt/m_predErr_sub_", linear_model, ".csv"))
m_predErr_sub_linear$X = NULL

m_predErr_sub_all = rbind(m_predErr_sub, m_predErr_sub_linear) 
m_predErr_all = rbind(m_predErr, m_predErr_linear)
m_predErr_all$WMSize = as.factor(m_predErr_all$WMSize)
levels(m_predErr_all$WMSize) =  c("low", "medium",  "high")
temp = m_predErr_all %>% filter(model == 'logarithmic') %>%summarise(abs_mmpredcv_err = abs(mmpredcv_err)) 

plt_Err_CV_all = ggplot(m_predErr_all, aes(abs(mmpredRP_err), abs(mmpredcv_err), group = interaction(model, Exp, WMSize), color = model, shape = Exp, size = WMSize)) + 
  geom_point() +
  geom_hline(yintercept = round(max(temp$abs_mmpredcv_err), 4)+0.0005, linetype='dashed')+
  xlab('Prediction error in the RP means (s)')+ ylab('Prediction error in CV')+colorSet3+
  scale_shape_manual(values = c(6, 7, 16, 17,8)) +
  theme_new+ 
  theme(legend.position = 'top')+
 labs(size = 'Memory Load')+
  guides(colour = guide_legend(order = 1, nrow=2,byrow=TRUE),
         shape = guide_legend(order =2, nrow=2,byrow=TRUE),
            size = guide_legend(order = 3, nrow=3,byrow=TRUE))

ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_Err_CV_all.png"), plt_Err_CV_all, width = 7, height = 4)
## Warning: Using size for a discrete variable is not advised.
plt_Err_CV_all
## Warning: Using size for a discrete variable is not advised.

m_predY_acc =  m_predErr_sub_all%>% 
  dplyr::group_by(Exp, model) %>% 
  dplyr::summarize(mmpredRP_rerr = mean(mpredRP_rerr)*100,
                   mmpredVar_rerr = mean(mpredVar_rerr)*100,
                   mmpredcv_rerr = mean(mpredcv_rerr)*100,
                   mmpredRP_acc = (1-mean(mpredRP_rerr))*100,
                   mmpredVar_acc = (1-mean(mpredVar_rerr))*100,
                   mmpredCV_acc = (1-mean(mpredcv_rerr))*100)
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
m_predY_acc